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SECTION  I 


* 


INTRODUCTION 

Owing  to  the  development  and  ever  greater  use  of 
advanced  composite  materials,  aerospace  structures  now 
evidence  the  most  rapid  evolution  since  the  introduction  of 
all  metal  aircraft.  High  strength  and  stiffness  characterize 
composites.  Such  efficiency  derives  from  the  intrinsic 
strength  of  the  constituent  materials  in  composites. 
However,  because  of  the  greater  complexity  of  their 
composition,  composites  demand  much  more  complex  structural 
evaluation  than  metals,  and  a  need  to  have  accurate  stress 
analysis  for  design  purpose  is  apparent. 

Delamination  is  a  critical  failure  mechanism  for 

laminated  composite  materials.  However,  before  delamination 

can  be  predicted,  it  is  essential  to  develop  analyses  to 
% 

accurately  estimate  the  interlaminar  stresses  that  cause 
delamination. 

Over  the  past  2  decades  considerable  effort  has  been 
devoted  to  finite  element  stress  analysis  of  structures  made 
of  composite  materials.  A  two-dimensional  analysis  using 
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plane  stress  or  plane  strain  assumptions  have  been  used  [1- 
5]  to  study  the  interlaminar  shear  strains  at  free  edge. 
Three-dimensional  finite  element  models  have  also  been 
developed  [6,7],  For  the  analysis  of  geometrically 
axisymmetric  composite  structures,  a  pseudo-three¬ 
dimensional  analysis  including  three  axisymmetric 
displacement  degrees  of  freedom  per  node  has  been  developed 
[8-10] .  In  this  approach,  the  complete  state  of  stress  and 
strain  is  modeled. 

Recently,  Chang  [11]  developed  a  more  sophisticated 
finite  element  model  based  on  continuous  displacement  as 
well  as  continuous  traction.  This  model  can  accurately 
predict  the  complete  stress  field  of  the  free  edge  stress 
problem  in  composite  laminate  coupons.  However,  it  is  quite 
expensive  in  computational  effort. 

The  Flight  Dynamics  Directorate  had  developed  a  model 
in  which  the  laminated  coupon  subjected  to  uniform  axial 
strain  in  longitudinal  direction  was  visualized  as  a  segment 
of  a  large  radius  axisymmetric  ring  under  uniform,,  radial 
pressure.  The  objective  of  the  present  research  was  to 
carefully  document  this  work,  improve  the  logic  of  the 
computer  code  to  make  it  more  efficient,  and  to  extend  its 
applicability  to  nonlinear  material  behavior. 


A  finite  element  computer  program  was  developed  for 
analysis  of  stresses  and  displacements  in  axisymmetric 
structures  made  of  laminated  composites.  Fiber  orientations 
at  an  angle  with  the  axis  were  allowed.  Three  degrees  of 
freedom  for  each  nodal  point  permit  modeling  of  the  complete 
state  of  stress  and  strain.  The  model  was  validated  by 
solving  representative  problems  of  cross-ply  and  angle-ply 
laminated  coupons  and  comparing  results  with  those  obtained 
using  the  continuous  traction  element  Q23  developed  by  Chang 
[11].  Section  II  contains  the  statement  of  the  axisymmetric 
problem  for  linearly  elastic  composite  materials.  In  Section 
III,  the  discretization  by  the  finite  element  method  is 
presented.  Section  IV  contains  application  of  the  computer 
program  developed  and  comparison  of  results  with  those 
obtained  using  Q23  finite  element  method.  The  procedure  was 
developed  specially  for  application  to  free-edge 
delamination  coupons.  However,  it  has  the  capability  to 
model  a  variety  of  structures  including  laminated  tubes, 
nozzles,  etc.  In  the  effort  described  in  this  report,  the 
material  was  assumed  to  be  linear  elastic.  Extension  of  the 
procedure  to  allow  for  nonlinear  material  behavior  is 
currently  in  hand  and  will  be  the  subject  of  another  report. 
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SECTION  II 


ELASTIC  AXISYMMETRICAL  PROBLEM 
FOR  CONTINUOUS  FIBER  COMPOSITE  MATERIALS 

The  equilibruim  equations,  strain  displacement 
relations  and  constitutive  relations  for  axisymmetric 
composite  materials  problem  are  listed  in  this  section. 
Derivation  of  the  stress-strain  matrix  in  material,  local 
and  global  coordinates  is  described. 

2.1  EQUILIBRIUM  EQUATIONS 

The  global  coordinate  system  for  axisymmetric  analysis 
is  a  cylindrical  system  where  r  is  the  radial  direction,  T 
is  the  circumferential  direction,  Z  is  the  axial  direction 
and  the  corresponding  displacement  are  u,  v,  w.  The 
equilibrium  equations  in  cylindrical  coordinates  are  [24]: 

2or/ar  +  (1/r)  a  ® rT/ a  T  +  *°rz/^z  +  (°r  “®T>/r  =  "fr 
3orz/<a r  +  (1/r)  <2oTz4)T  +  Joz/dz  +  orz/r  =  -fz 

a<>rT/^r  +  (Vr)  a«T/aT  +  ^°Tz/«JZ  +  2  arT'/r  =  -fT 

(2.1) 


where  fr,  fz,  fT  are  components  of  body  forces  in  r,  z,  T 
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directions  respectively.  For  axi symmetric  problems  the 
stresses  are  independent  of  the  tangential  coordinate  T. 
Therefore,  the  equilibrium  equations  reduce  to: 

^or/^r  +  boTZ/} z  +(or  -oT)/r  =  -fr 

^or2/^r  +  *oz/d z  +  orz/r  =  -fz  (2.2) 

^orT/3  r  +  ^oTz/^z  +  2  orT/r  =  -fT 

2.2  STRAIN-DISPLACEMENT  RELATIONS 

The  general  strain-displacement  relations  in 

cylindrical  coordinates  are  given  by  [24]: 

er  =  bu/)r  eT  =  u/r  +  (l/r) ft  v/^T) 

ez  =  dw/^z 

2erT  =  (l/r)(*u/aT)  +  (^v/ar)  -  v/r  (2.3) 

2erz=^u/az  +iw/Jr 

2ezT  =  ^v/*z  +  (l/r)  (^ w/^T) 

For  axisymmetric  analysis,  the  three  components  of 
displacement  are  all  independent  of  T.  Thus,  the 

relationships  given  by  (2.3)  reduce  to: 

er  =  bn/br  ,  eT  =  u/r 

ez  =  fcw/dz  ,  2erz  =  ^u/^z  +  ^w/^r  (2.4) 

2czT  =  bv/bz  ,  2erT  =  -  v/r 
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2.3  STRESS-STRAIN  RELATIONS  FOR  ORTHOTROPIC  COMPOSITE 


MATERIALS 

Composite  structures  made  of  oriented  continuous  fibers 
are  usually  constructed  so  that  the  high-strength,  high- 
stiffness  fibers  are  embedded  in  a  weaker  matrix  material. 
Though  the  material  is  heterogeneous,  it  is  often 
characterized  as  an  equivalent  homogeneous  material  for  use 
in  analysis  [10],  [12].  In  the  present  analysis  the  material 
is  assumed  to  be  linear,  homogeneous  and  orthotropic. 
Because  of  the  material  axis  not  being  parallel  to  the 
global  reference  frame  (Figure(l)),  the  stress-strain 
relations  may  appear  to  be  anisotropic  when  written  with 
reference  to  the  global  axes. 

2.3.1  Stress-Strain  Relationship  in  Material  Coordinates 

In  the  plane  of  fibers  the  material  is  orthotropic, 
and,  therefore,  in  three-dimensional  configuration,  there 
are  nine  independent  elastic  constants.  Let  (1,2,3)  be  the 
material  directions.  Using  'reduced'  vector  representation 
for  the  symmetric  stress  and  strain  tensors,  the  stress- 
strain  relationship  can  be  written  in  matrix  form  as  [12]: 

<°>n,  -  [C]m  {6}ffl 

Explicitly: 
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(2.5) 


Where : 

Clls=  E11  *  1  “  ^23^32  ^  /D 


C22  E22(1 

-  *i31»i13)  /° 

C33~  E33  ^ 1 

^12^21)/° 

C44=  G23 

C55=  G13 

C66  G12 

(2.6) 

C12=  ^21+ll31^23^  Ell//D  =  ^12+>113^32)E22/D 

C13=  (p31+ll21M32)  Ell/D  =  ^13+^12ll23)E33/D 

C23=  (u32+J112M31)E22/,D  =  ( JA23  +  >i21lA13 J  E33/ D 

D  ~  1  ^12^21”  ^23^32”  ^31^13"  2>121,132,113 

and  E^=  Modulus  of  elasticity  in  the  i  direction  i=l,2,3 

G^j=  Shear  modulus  in  the  i-j  plane,  i,j=l,2,3 

Poisson's  ratio  in  the  i-j  plane  i,j=l,2,3 
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2.3.2  Stress-Strain  Relationship  in  Local  Coordinates 


In  order  to  combine  the  standard  behavior  of  various 
laminae  to  get  the  properties  of  the  angle-ply  laminate,  it 
is  necessay  to  rewrite  the  constitutive  relations  in  a 
uniform  system  of  coordinates.  To  do  this  one  has  to 
establish  the  relationships  between  stresses  and  strains  in 
the  material  coordinate  system  with  corresponding  quantities 
in  the  local  coordinate  system. 

The  analysis  presented  herein  is  for  an  angle  of 
rotation  a  of  the  fiber  direction  (1)  as  measured  from  the 
Z  (axial)  direction.  Thus  if  B  is  the  fiber  orientation 
measured  from  the  tangential  direction, 
a  -  B-  90.  (Figure  (2)).  Let: 


a  1 

H 

_ 1 

°z 

“ 

*z 

°2 

°T 

e2 

®T 

°3 

=  [J] 

°rz 

f 

*3 

=  [T] 

er 

°32 

°rT 

2e32 

2erT 

°31 

°rz 

2e31 

2erz 

^12_ 

_ °zT_ 

J:ei2_ 

_^®ZT_ 
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or  symbolically 


(°>m  =  [J]  {°}L  ;  {e}m  =  [T]  { e }L 

where  [T] ,  [J]  are,  respectively,  the  strain  and  the  stress 
transformation  matrix  due  to  a  rotation  of  the  axis  through 
an  angle  a  and  subscripts  m  and  L  refer  to  material  and 
local  coordinate  systems  respectively.  Further,  since  strain 
energy  is  invariant  due  to  rotation  of  axis  we  have: 

[J]"1  =  [T]t 

Let  {t}  be  the  traction  vector,  and  {n}  the  normal  at 
a  certain  point  in  the  material  ,  then  by  Cauchy's  stress 


principle,  (t)  =  [o] 

{n} , 

where  [o] 

is  the 

stress  matrix: 

— 

°11 

°12 

°13  _ 

[o]  = 

°12 

°2  2 

°23 

(2.8) 

”31 

°32 

°33  _ 

Due  to  a  rotation 

of 

axis  about  r. 

the  new  traction  and 

normal  vector  are  given  as: 


{t>m  =  [L]  {t}L  { n)m  »  [L]  { n)L  (2.9) 


Where  [L]  is  the  transformation  matrix  due  to  a  rotation  of 
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axis  Z,  T  about  r  by  an  angle  «.  Explicitly: 


in  which  m  =  Cos  «,  n  =  Sin  a. 


(2.10) 


The  new  stress  matrix  [o]^  in  the  local  coordinates  is 
then  given  by: 

[«3i  =  [L]"1  («]  [L], 

Comparing  term  by  term  and  using  the  'reduced*  (vector) 
representation  for  stress,  we  get  the  stress  transformation 
matrix  relating  stress  in  local  coordinates  to  stress  in 
material  coordinates  [18]. 


(2.11) 
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The  matrix  [J]  can  be  easily  computed  from  [J]-1.  It 
corresponds  to  the  transformation  due  to  a  rotation  -«  about 
r  axis.  Therefore,  [J]  can  be  found  by  replacing  *n'  by  ' - 
n'  in  [J]-1  i.e., 


(2.12) 

The  stress-strain  relationship  in  local  coordinates  is 
then  determined  as  follows: 


We  have. 

t°>m 

<«>m 

(o}L 

1! 

r— i 

<°>m' 

(2 

.13) 

<*>m 

=  [T] 

(e>L 

then 

(°>L  : 

=  [T]T 

CC]B  [T] 

{e>L 

(2 

•  14) 

Therefore, 

the 

stress 

-strain 

relationship 

in  local 

coordinates 

is  given  by: 

[C]L  = 

[T]T 

[Clm  CT] 

Explicitly: 
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For  the  case  of  a  ring  in  which  the  local  and  global 
axes  do  not  coincide,  it  is  necessary  to  transform  the 
stress-strain  relations  derived  in  local  system  of  axes  to 
global  coordinates.  We  note,  however,  that  for  cylindrical 
objects  global  and  local  axes,  in  general,  coincide.  The 
matrix  relating  stresses  to  strains  after  rotation  of  axis 
about  the  tangential  axis  is  derived  following  the  same 
approach  used  in  the  preceding  section.  Let: 

{<7>l  -  [Jx]  (o)a 

(«>i  -  [TJ  { c )c  (2.17) 

(t}L  =  [L]  {t)c 
(n}L  =  [L]  {n}c 

where  the  subscripts  G,  L  denote  the  global  and  the  local 
coordinates  systems  respectively,  and 


(2.18) 


in  which  M=cos  a',  N=sin  a' ,  and  a'  is  the  angle  of  rotation 
of  r-Z  plane  about  the  tangential  T  axis.  Then  the  inverse 


[L]  = 


M 

0 

N 


0  -N 

1  0 

0  M 
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of  the  stress  transformation  matrix  relating  stress  in 
global  coordinate  to  stress  in  local  coordinates  is: 


M2 

0 

N2 

0 

2MN 

0 

0 

1 

0 

0 

0 

0 

[Ji]"1-  [Tx] 

T_ 

N2 

0 

M 

0 

-2MN 

0 

0 

0 

0 

M 

0 

-N 

-MN  0 

MN 

0 

m2-n2 

0 

_  0 

0 

0 

N 

0 

M 

The  stress-strain 

matrix 

in  the 

global 

coordinate  : 

computed 

as: 

[C]G  = 

[tx]t 

[C]L  [T^. 

Explicitly: 

— 

DC11 

DC12 

DC13 

DC14 

DC15 

DC16 

DC12 

DC22 

DC23 

DC24 

DC25 

DC26 

[C]G  - 

DC13 

DC23 

DC33 

DC34 

DC35 

DC36 

DC41 

DC42 

DC43 

DC44 

DC45 

DC4  6 

t 

DC15 

DC2  5 

°C35 

DC45 

DC55 

DC56 

DC16 

DC2  6 

°C36 

DC46 

DC56 

DC66 

(2.19) 


2.20) 
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where 


DC11  =  CC11  M4  +  2  CC13  M2N2  +  CC33  N4  +4M2N2  CC55 

DC12  =  CC12  M2  +  CC23  N2 

DC13  =  (CC1X  +  CC33  -  4CC55)  M2N2  +  CC13  (M4  +  N4) 

DC14  =  2  M2N2  CC45  -M2N  CC16  -N3  CC36 

dc15  =-mn3  cc11+m3n  CC33+MN (m2-n2)  cc13+(m2-n2)2Mn  cc55 

DC16  =  MN2  (2  CC45  +CC16)  +  M2N  CC36 
DC22  =  CC22 

DC23  =  N2  CC12  +  M2  CC23 
DC24  =  -N  CC26 

DC25  =  CC12  +  ^  CC23 

DC26  =  M  CC26 

DC33  =  CC33  M4  +  2  CC13  M2N2  +  CCX1  N4  +4M2N2  CC55 
DC34  =  2  M2N  CC45  -M2N  CC26  -N3  CC16 

DC35  =— MN3  CC11+M3N  CC33-  MN(M2-N2)  CC33- (M2-N2) 2MN  CC55 
DC36  =  MN2  (~2  CC45  +CC26)  +  M3  CC26 
DC44  =  M2  CC44  +  N2  CC66 

DC45  =  MN2  (CC16  -  CC36)  +  (M2  -  N2)M  CC45 
DC46  =  MN  CC44  -  MN  CC66 

DC55  =  (M2N2)(CC11  -2CC13  +  CC33)  +  (M2  -  N2)2  CC55 
DC56  =  (M2N  "  N3)  CC45  +  M2N  (~CC16  +  CC36) 

DC66  =  N2  CC44  +  M2  CC66 

(2.21) 


17 


SECTION  III 


DISCRETIZATION  BY  THE  FINITE  ELEMENT  METHOD 

3.1  BASIC  CONCEPTS: 

The  finite  element  method  is  a  general  technique  for 
construction  of  approximate  solutions  to  boundary  value 
problems.  The  method  essentially  involves  dividing  the 
spatial  domain  of  interest,  R,  into  a  finite  number  of 
nontrivial  discrete  elements,  and  using  variational  concepts 
to  construct  an  approximation  of  the  solution  over  the 
collection  of  finite  elements. 

Therefore,  the  principal  ingredients  of  the  finite 
element  method  for  constructing  approximate  solutions  of 
problems  are: 

1.  The  construction  of  a  finite  element  mesh  and 
selection  of  a  set  of  piecewise-polynomial  basis  functions 
defined  on  the  mesh  satisfying  the  required  continuity 
conditions.  The  elements  are  assumed  to  be  interconnected  at 
a  discrete  number  of  nodal  points  situated  on  their 
boundaries. 

2.  The  construction  of  an  approximation  of  the 
variational  formulation  of  the  problem,  setting  up  of 
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element  properties,  assemblage  of  these  into  a  system  of 
algebraic  equations  with  nodal  values  of  the  dependent 
variables  or  their  derivatives  up  to  a  certain  order  as  the 
unknown  quantities. 

3.  Solution  of  the  algebraic  system. 

3.2  FORMULATION 

Sandhu  and  Pister  [15,16]  and  Sandhu  and  Salaam  [17] 
used  a  generalization  of  Mikhlin's  basic  variational  theorem 
to  extend  its  application  to  linear  coupled  initial  boundary 
value  problem.  Specialization  of  an  "extended"  variational 
principle  leads  to  the  well-known  potential  energy  theory 
which  has  been  used  in  developing  the  finite  element  model 
described  in  this  report.  For  elastostatics  problem, 
assuming  that  the  strains  exactly  satisfy  the  strain- 
displacement  relationship,  and  the  displacement  field  is 
continuous  over  the  spatial  region  of  interest  R  and 
satisfies  the  prescribed  displacement  boundary  conditions  on 
portion  S^^  of  the  boundary  d  R  of  R,  the  potential  energy  is: 

I(u>  =  <“i,j'  Eijkl  uk,l>R  "  2<ui'  fi>R 

—  2<u^,t^>g2  (3.1) 

Here  u^,  f^,  t^  are  components,  respectively  of  the 

displacement  vector,  the  isothermal  elasticity  tensor,  body 
force  vector,  and  the  specified  tractions  on  boundary  S2. 
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<  ,  >R,  <  ,  >S2  denotes  inner  products  of  pairs  of  functions 
evaluated  over  the  domain  indicated  by  the  subscript.  S1,  S2 
are  complementary  subsets  of  ^  R. 

Vanishing  of  the  variation  of  I(u)  with  respect  to  the 
displacement  components  u^  implies  the  satisfaction  of 
equilibrium  equations  over  R  and  on  S2.  The  finite  element 
approximation  consists  of  applying  the  theorem  to  the 
approximate  functional  evaluated  by  using  the  approximate 
values  of  the  functions  defined  by  the  finite  element 
interpolation  scheme. 

Taking  the  variation  of  the  discretized  form  of  the 
functional  yields  the  system  of  algebraic  equations: 

[K] { r }  =  {F} 
where : 

[K] :  The  system  stiffness  matrix. 

{r}:  Nodal  point  displacement  vector. 

{F}:  Nodal  point  load  vector. 
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3.3  INTERPOLATION  FUNCTIONS  FOR  THE  MODIFIED  BILINEAR 

ISOPARAMETRIC  QUADRILATERAL  ELEMENT  OM4 

3.3.1  Bilinear  Interpolation  Function  for  Isoparametric 
Quadrilateral  Element 

The  element  used  in  the  analysis  is  an  arbitrary 
quadrilateral.  The  displacements  are  bilinearly  interpolated 
over  each  element.  The  interpolation  functions  are  expressed 
as: 

Nt  =  1/4  (1-s) (1-t) 

N2  -  1/4  (1+s) (1-t)  (3.2) 

N3  =  1/4  (1+S) (1+t) 

N4  *  1/4  (1-s) (1+t) 

where  s,t  is  a  set  of  coordinates  such  that  the  corners  of 
the  quadrilateral  are  (-1,-1),  (-1,1),  (1,1),  (1,-1). 

3.3.2  Modified  Bilinear  Isoparametric  Element 
Selective  reduced  integration  was  introduced  by  Taylor 

et  al.,  [19]  for  evaluation  of  the  stiffness  matrix  for 
bending  problems.  Salaam  [18]  reported  remarkable 
improvement  in  numerical  results  for  some  problems  where 
selective  reduced  integration  had  been  used.  This  was 
physically  explained  [19]  as  resulting  from  the  excessive 
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transverse  shear  strain  energy  contribution  when  using 
isoparametric  elements.  As  a  remedy  for  this  problem,  a 
constant  shear  strain  was  imposed  on  the  element  by 
assigning  the  shear  strain  associated  with  the  center  of  the 
element  to  the  entire  element.  The  modified  (constant  shear) 
element  so  obtained  was  labelled  as  QM4  and  was  used  in  the 
finite  element  code  developed  in  this  research  program. 
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SECTION  IV 


APPLICATION  TO  DELAMINATION  SPECIMENS 

4.1  INTRODUCTION 

Although  the  computer  program  is  capable  of  obtaining 
approximate  stresses  and  deformations  in  any  axisymmetrical 
laminated  solid,  the  primary  motivation  of  the  development 
was  to  obtain  an  approximation  of  stress  fields  in  a  free- 
edge  delamination  coupon.  The  program  was  validated  through 
application  to  a  four-ply  symmetric  laminate.  Results  were 
compared  to  numerical  solutions  obtained  using  a  continuous 
traction  finite  element  formulation  [11].  The  computer 
program  was  also  used  to  solve  for  stresses  and  deformations 
in  a  22-layer  symmetric  angle-ply  laminate. 

4.2  MATHEMATICAL  MODEL 

The  straight  coupon  subjected  to  a  uniform  inplane 
strain  in  the  longitudinal  direction  is  viewed  as  a  short 
segment  of  a  large  radius  ring  subjected  to  a  radial 
pressure.  Figure  (3)  shows  the  correspondence  between  a  ring 
segment  and  a  short  straight  coupon.  For  a  sufficiently 
large  radius  of  the  ring,  a  uniform  inner  pressure  will 
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R,  Z,  T  in  a) 
corresponds  to 
Z,  Y,  X  in  b) 


/  (b) 

X 

/ 

Figure  3:  Correspondence  Between  Ring  Specimen  a) 
and  Strip  Specimen  b) 


yield  a  reasonably  uniform  tangential  strain  over  all 
elements  of  a  section.  Also,  for  large  radius  the  magnitude 
of  the  radial  internal  pressure  will  be  small  in  comparison 
with  the  other  stresses  in  the  ring  and,  therefore,  will  not 
be  of  much  influence  in  determining  failure. 

To  insure  accuracy  of  results,  different  values  of  the 
inner  radius  were  tried.  It  was  found  that  for  a  radius 
equal  to  104  times  the  thickness  of  the  specimen,  the 
tangential  strain  was  uniform,  up  to  a  four-digit  accuracy, 
over  all  elements  of  the  section.  Figure  (4)  shows  through¬ 
thickness  distribution  of  T-strain  near  the  free-edge  for 
radii  104,  103  and  102  times  the  thickness. 

In  general,  axisymmetric  loading  of  a  composite  ring 
will  cause  circumferential  displacement  due  to  the  composite 
material  behavior.  Thus,  when  the  hoop  is  loaded 
axisymmetrically  with  internal  pressure,  a  circumferential 
strain  is  generated  in  addition  to  the  radial  deformation 
and  shear  strain  over  planes  parallel  t  the  axis  [10], 
Therefore,  to  allow  the  complete  state  of  stress  and  strain 
to  be  modeled,  three  displacement  degrees  of  freedom, 
independent  of  the  circumferential  coordinate,  need  to  be 
associated  with  each  node:  u,  axial;  v,  circumferential;  w, 
radial.  Two  computer  programs  were,  therefore,  developed. 
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One  allcr  2d  for  all  the  three  displacement  degrees  of 
freedom  and  the  other  allowed  only  the  radial  and  axial 
displacements.  In  the  two-dimensional  analysis,  rT  and  zT- 
strains  are  assumed  zero  while  the  corresponding  stresses 
can  have  nonzero  values,  depending  upon  the  material 
properties  and  fiber  directions. 

4.3  FOUR-PLY  LAMINATES 

4.3.1  Geometry.  Material  Properties  and  Finite  Element 
Mesh 

Ring  specimens  with  fiber  orientation  [0/90]s  and  [45/- 
45] s  under  uniform  internal  pressure  were  considered  in  the 
present  analysis.  The  relation  between  laminate  width  and 
thickness  was  2b  =  16h,  as  in  references  [11],  [23].  Each 

ply  was  idealized  as  a  homogeneous,  elastic  orthotropic 
material  with  the  following  properties  [11]: 

Elx  =  20  *  106  psi 

E22  =  E33  =  2.1  *  10®  psi 

G32  “  G^3  ”  G2 3  =  0.85  *  10®  ps 1 

^12  =  ^13  =  ^23  =  0 ' 21 

The  subscripts  1,2,3  correspond  to  material  coordinates 
in  the  longitudinal,  transverse  and  thickness  directions 
respectively.  Each  physical  layer  of  thickness  h  was 
subdivided  into  sublayers,  the  total  number  of  sublayers  is 


denoted  by  N.  The  mesh  used  was  refined  near  the  free  edge 
in  order  to  obtain  more  reliable  results  in  regions  where 
steeply  varying  stresses  are  expected.  A  288  (N=16)  element 
model,  as  shown  in  Figure  (5),  was  used  to  discretize  a 
typical  half  section  of  the  ring.  This  corresponds  to  a  144- 
element  model  for  a  quadrant  as  used  with  continuous 
traction  Q23  element  in  [11].  The  r,  z,  T  axes  in  the  model 
correspond  to  the  z,  y,  x  axes  in  the  cartesian 
representation  of  the  coupon.  The  edge  z=b  is  the  free-edge 
and  R=0,  R=4h  respectively  are  the  bottom  and  top  surface  of 
the  specimen.  The  value  R^,  shown  in  Figure  (5)  ,  denotes  the 
distance  from  the  central  plane  of  section  to  the  center  of 
the  ith  sublayer.  R-^  is  the  closest  to  the  midplane  of  the 
specimen  and  Rg  is  closest  to  the  top  surface. 

The  internal  pressure  applied  was  1108  (psi)  and  290.53 
(psi)  for  the  cross-ply  and  the  angle-ply  laminate, 
respectively.  This  pressure  yields  a  unit  average  axial 
strain,  corresponding  to  the  uniform  axial  strain  applied  in 
the  continuous  traction  Q23  analysis  [11]. 
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Radius  =  104  *  Thickness 
Figure  5:  288  Element  Mesh  (16x18) 
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4.3.2  Cross-Plv  Laminates 

The  cross-ply  laminate  was  analyzed  using  both  the  two 
degree-of-freedom  computer  program  as  well  as  the  pseudo- 
three-dimensional  approach  and  the  results  compared.  For 
cross-ply  laminates,  the  transverse  Tr  and  Tz  stresses  were 
negligible  compared  to  the  other  stresses  in  the  pseudo- 
three-dimensional  analysis  and  were  nearly  zero  in  the  two- 
dimensional  analysis.  Further,  the  four  stress  components 
or,  o2,  oT,  or2  were  only  slightly  affected  by  the 
introduction  of  the  third  degree  of  freedom.  Figures  (6) , 
(7)  show  the  distribution  of  R-stress  along  R1  and  Rz-stress 
along  R4  for  the  cross-ply  laminate.  The  figures  indicate 
excellent  agreement  between  results  obtained  using  three 
degrees  of  freedom  for  each  node  and  those  obtained  using 
two  degrees  of  freedom  for  each  node.  These  sections  were 
selected  for  the  comparison  as  representatives.  Similar 
agreement  was  also  obtained  for  stress  distributions  at 
other  sections. 

Because  of  the  symmetry  in  the  laminate  about  the  Z=0 
axis,  only  half  the  section  was  considered.  Numerical 
results  based  on  two-dimensional  axisymmetric  analysis  using 
QM4  element  were  compared  with  results  based  on  the 
continuous  traction  Q23  element  described  in  [11].  The  same 
mesh  size  was  used  in  both  the  analyses.  The  numerical 
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Figure  6:  Distribution  of  R-Stress  Along  R1  of 

[0/90]s  Laminate:  2-D  and  Pseudo  3-D  Analysis 
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Figure  7:  Distribution  of  Rz-Stress  Along  R4  of 

[ 0/90 ] 8  Laminate:  2-D  and  Pseudo  3-D  Analysis 
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result  using  Q23  [11]  has  shown  good  agreement  with  Pagano's 
analytical  solution  [23]  which  is  based  on  a  generalization 
of  Reissner's  theory.  Figures  (8)  through  (23)  compare 
stress  fields  at  specific  locations  for  the  cross-ply 

♦ 

laminates  using  Q23  and  the  axisymmetric  analysis. 

v. 

Distribution  of  R-stress  along  the  center  of  the  first 
sublayer  beyond  the  central  plane  (at  Rl)  of  the  [0/90]s  is 
shown  in  Figure  (8)  .  Solutions  obtained  from  the 
axisymmetric  analysis  practically  coincide  with  Q23  element 
solution  over  the  entire  width  of  the  laminate.  The 
distribution  of  stresses  indicates  a  sharp  rise  near  the 
free  edge  which  agrees  with  Pagano's  solution  [11], [23]. 
Figures  (9),  (10)  show  the  distribution  of  R-stress  at 

center  of  the  last  sublayer  of  0°  ply  (at  R4)  and  the  first 
sublayer  of  the  90°  ply  (at  R5)  respectively.  At  the 
interface  of  the  0  and  90  degree  plies,  a  stress  singularity 
is  expected  at  the  free  edge  due  to  the  discontinuity  in 
elastic  properties.  The  stress  distributions,  however,  are 
only  obtained  at  the  centers  of  adjacent  elements  and  can 
only  be  compared  there.  It  is  noted  that  the  stress  value 
near  the  free-edge  obtained  using  axisymmetric  analysis 

% 

overestimated  the  stress  at  R4  and  underestimated  it  at  R5 
as  compared  with  Q23  solution. 
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Figure  8:  Distribution  of  R-Stress  Along  Rl 
of  [ 0/90 ] s  Laminate. 
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Figure  10:  Distribution  of  R-Stress  Along  R5 
of  [0/90]s  Laminate. 
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A  discrepancy  in  the  distribution  of  R-stress,  for  the  two 
procedures,  at  the  center  of  the  last  (nearest  to  the 
surface)  sublayer  (R8)  and  the  sublayer  near  to  it  is 
observed  in  Figures  (11)  and  (12).  This  could  possibly  be 
because,  in  the  present  procedure,  the  traction-free 
boundary  condition  is  not  enforced.  The  difference  in 
stresses  is  less  severe  away  from  the  free  surface  as  shown 
in  the  distribution  of  R-stress  at  R7  (Figure  (11))  than 
that  at  R8  (Figure  (12)). 

Figure  (13)  shows  through-the-thickness  distribution  of 
R-stress  near  the  free  edge  (at  Z/B  =  0.99)  of  the  laminate. 
In  the  vicinity  of  the  interface,  a  difference  in  the 
solution  of  continuous-traction  element  and  the  axisymmetric 
analysis  exists.  It  was  felt  that  this  could  be  due  to  the 
finite  element  mesh  not  being  fine  enough  to  approximate  the 
steeply  varying  stresses  associated  with  abrupt  change  of 
material  properties.  To  explore  this  further,  the  effect  of 
mesh  refinement  was  studied.  This  is  described  later  in  this 
section. 

Values  of  Rz-stress  along  the  center  of  the  first 
sublayer  beyond  the  center  plane  of  the  laminate  (along  Rl)  , 
along  R4  and  along  R5  are  shown  in  Figures  (14,15,16) 
respectively.  The  result  showed  satisfactory  agreement  with 
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those  obtained  using  Q23  element.  However,  a  discrepancy  is 
observed  near  the  free  edge  in  Figure  (16)  .  This  is 
apparently  because  the  continuous  traction  Q23  element 
exactly  satisfies  the  traction-free  boundary  condition, 
while  axisymmetric  displacement  formulation  does  not.  A  more 
refined  mesh  near  the  free  edge  is  likely  to  give  more 
satisfactory  result.  This  was  investigated  and  is  described 
later  in  this  section. 

Figures  (17)  and  (18)  show  the  distribution  of  Rz- 
stress  along  R7  and  along  the  center  of  the  last  sublayer 
(at  R8) .  Contrary  to  the  distribution  of  R-stress  at  these 
locations,  the  RZ-stress  distributions  are  in  good  agreement 
with  those  obtained  from  Q23  element. 

Through-the-thickness  distribution  of  Rz-stress  near 
free  edge  (at  Z/b  =0.99)  is  shown  in  Figure  (19)  .  The 
results  show  significant  difference  with  those  for  the  Q23 
analysis  specially  near  the  interface  of  0°/90°  layers. 

Distribution  of  Z-stress  along  Rl,  R4  and  R8  are  shown 
in  Figures  (20,  21,  22)  respectively.  The  results  along  Rl 
and  R8  obtained  using  the  axisymmetric  analysis  agreed  well 
with  the  Q23  element  solution  across  the  entire  width  of  the 
laminate.  However,  a  difference  in  stress  distribution  along 
R4  is  observed  near  the  free  edge  where  a  steeply  varying 
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Figure  14:  Distribution  of  RZ-Stress  Along  R1 
of  [0/90] s  Laminate. 
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Figure  15:  Distribution  of  RZ-Stress  Along  R4 
of  [0/90]s  Laminate. 
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Figure  16:  Distribution  of  RZ-Stress  Along  R5 
of  [ 0/90 ] s  Laminate. 
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Figure  17:  Distribution  of  RZ-Stress  Along  R7 
of  [ 0/90 ] s  Laminate. 
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Figure  18:  Distribution  of  RZ-Stress  Along  R8 
of  [ 0/90 ] s  Laminate. 
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stress  exists.  Figure  (23)  shows  the  distribution  of  T- 
stress  along  R8.  The  figure  indicates  satisfactory  agreement 
between  the  axisymmetric  analysis  and  the  Q23  element 
solution. 

4.3.3  Anale-Plv  Laminate 

A  ring  specimen  with  fiber  orientation  [45/-45]s  under 
uniform  internal  pressure  was  analyzed.  Using  the  two 
axisymmetric  analysis  approaches,  it  was  found  that  the 
introduction  of  the  circumferential  degree  of  freedom  was 
necessary  for  modeling  the  complete  state  of  strain  and 
stress.  The  distribution  of  Z-stress  along  R1  and  Rz-stress 
along  R4,  Figures  (24)  and  (25)  respectively,  indicate  that 
the  two-dimensional  analysis  gave  negligible  stresses  as 
compared  to  results  obtained  using  the  pseudo  three- 
dimensional  analysis.  The  Tr-stress  had  values  of  the  order 
106  psi  in  the  pseudo  three-dimensional  analysis,  comparable 
to  the  z,  rz  and  r-stress  components,  while  they  were  zero 
in  the  two-dimensional  axisymmetric  analysis.  Figure  (26) 
shows  the  distribution  of  Tz-stress  along  R1  and  indicates 
that  the  two-dimensional  analysis  gave  a  practically 
constant  variation  of  stress  and  did  not  satisfy  the  stress- 
free  boundary  condition  at  the  free  edge  while  the  pseudo- 
three  dimensional  analysis  did. 
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Figure  24:  Distribution  of  Z-Stress  Along  R1  of 

[45/-45]s  Laminate:  2-D  and  Pseudo  3-D  Analysis 
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Figure  25:  Distribution  of  Rz-Stress  Along  R4  of 

[ 4 5/-4 5 ] s  Laminate:  2-D  and  Pseudo  3-D  Analys 


For  angle-ply  laminate,  in  contrast  to  the  case  of  the 
cross-ply  laminates  with  material  axes  parallel  to  the 
tangential  plane,  the  symmetry  about  the  Z  =  0  plane  does 
not  exist.  Both  half-section  and  full-section  of  the  ring 
specimen  were  analyzed  using  the  same  mesh  and  results  were 
compared.  Figures  (27)  and  (28)  show  the  distribution  of  Tz- 
stress  and  Z-stress  along  R1  for  full  and  half  section 
analysis.  Even  though  the  results  agree  near  the  free 
surface,  the  stresses  near  midsurface  were  quite  different. 
Therefore,  for  analysis  of  angle-ply  laminates, 
discretization  of  the  full  section  is  necessary. 

Results  for  the  pseudo-three-dimensional  analysis  were 
compared  with  the  numerical  solution  obtained  using  Q23 
element  [11].  The  same  mesh  size  was  used  in  both 
procedures.  Figures  (29)  through  (44)  compare  stress  fields 
at  specific  locations  for  the  [45/-45]s  laminate  using  the 
Q23  and  the  present  axisymmetric  analysis.  The  mesh  used  was 
refined  near  the  free  edge.  Each  physical  layer  was 
subdivided  into  sublayers,  the  total  number  of  sublayers  is 
denoted  by  N.  Full  section  was  used  for  the  axisymmetric 
analysis,  corresponding  to  288  (16x18),  (N=16)  mesh  for  half 
the  section  of  the  ring  (i.e.,  the  total  number  of  elements 
used  was  576)  .  Figures  (29)  ,  (30)  and  (31)  show  the 
distribution  of  Tz-stress  along  Rl,  R4  and  R8  respectively. 
The  results  along  Rl  and  R8  obtained  using  axisymmetric 
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Figure  27:  Distribution  of  Tz-Stress  Along  R1  of 
[45/-45]s:  Full  and  Half  Section  Analysis 
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Figure  31:  Distribution  of  Tz-Stress  Along  R8 
of  [ 4 5/-4 5 ] s  Laminate 


analysis  agreed  well  with  Q23  element  solution  across  the 
entire  width  of  the  laminate.  The  distribution  of  other 
stress  components  along  R1  were  also  in  satisfactory 
agreement  with  Q23  element  solution.  A  difference  in  the 
stress  distribution  along  R4  is  observed  near  the  free  edge 
where  a  steeply  varying  stress  exists. 

A  comparison  of  the  TR-shear  stress  distribution  along 
R4  (Figure  (32))  indicated  that  the  axisymmetric  analysis 
solution  had  sharp  rise  toward  the  free  edge  similar  to  Q23 
solution.  Satisfactory  agreement  was  observed  between  these 
two  solutions  for  stress  across  the  width  except  at  the 
free-edge  boundary  where  axisymmetric  analysis  overestimated 
the  stress  values  as  compared  to  Q23.  Figure  (33)  shows  the 
distribution  of  TR-stress  along  R8.  The  results  agreed  well 
except  near  the  free  edge.  This  could  be  due  to  the  finite 
element  mesh  not  being  fine  enough  to  approximate  the 
steeply  varying  stress  near  the  free  corner.  Figure  (34) 
shows  through-thickness  distribution  of  TR-stress  near  the 
free  edge  (at  Z/B=0.99).  Singular  behavior  of  TR-stress 
which  is  highly  localized  near  the  interface  of  45/-45 
layers  is  noticeable.  The  axisymmetric  analysis  gave  higher 
value  of  stresses  than  the  Q23  element  in  the  vicinity  of 
the  interface.  Distribution  of  R-stress  along  R4  is  shown  in 
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Figure  32:  Distribution  of  TR-Stress  Along  R4 
of  [ 45/-45]s  Laminate 
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Figure  34:  Through-Thickness  Distribution  of  TR- 
Stress  Near  Free  edge  (Z/B=0.99)  of 
[45/-45]s  Laminate.  Q23  and  axisyiranetric 
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Figure  (35) .  The  results  showed  reasonable  agreement  with 
those  obtained  using  Q23  element.  Through-thickness 
distribution  of  R-stress  near  the  free  edge  (Figure  (36)) 
shows  a  steeply  varying  stress  distribution  in  the  vicinity 
of  the  interface  of  45/-45  layers. 

Figure  (37)  shows  the  distribution  of  Rz-stress  along 
R4 .  A  very  close  agreement  was  generally  observed  between 
solutions  except  at  the  last  value  near  the  free  edge  where 
a  steeply  varying  stress  exists.  Figure  (38)  shows  through¬ 
thickness  distribution  of  Rz-stress  at  Z/B=0.99.  A  steeply 
varying  stress  distribution  associated  with  abrupt  change  of 
material  properties  is  also  observed  in  the  vicinity  of  the 
interface  of  45/-45  layers.  However,  the  axisymmetric 
analysis  gave  a  lower,  as  compared  to  Q23  analysis,  estimate 
of  the  stresses  near  the  interface.  A  discrepancy  in  the 
distribution  of  Z-stress  along  R4  and  R8,  between  the  two 
procedures,  near  the  free  edge  is  observed  in  Figures  (39) 
and  (40) . 

Distribution  of  T-stress  along  R8  is  shown  in  Figure 
(41) .  Contrary  to  the  T-strain  distribution,  the  T-stress 
distribution  is  not  uniform  along  a  sublayer.  This  is 
because  it  depends  on  r,  z,  Tz  and  T-strains.  Figures  (42- 
44)  show  the  distributions  of  r,  z,  and  Tz-strains  near  the 
free  surface  of  the  coupon. 
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Figure  36:  Through-Thickness  Distribution  of 
R-Stress  Near  Free  edge  (Z/B=0.99) 
of  [45/-45]s  Laminate 
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Figure  37:  Distribution  of  Rz-Stress  Along  R4 
of  [45/-45]s  Laminate 
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Figure  38:  Through-Thickness  Distribution  of  Rz- 
Stress  Near  Free  edge  (Z/B=0.99) 
of  [ 4 5/-4 5 ] s  Laminate 
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4.3.4  Effect.  Of  Mesh  Refinement 


Analysis  of  the  four-ply  laminate  was  done  in  this 
study  using  a  288-element  model  for  half  section  analysis, 
locally  refined  near  the  free  edge  in  the  Z-direction.  To 
study  the  effect  of  local  mesh  refinement  near  the  free  edge 
on  the  solutions,  comparisons  were  made  for  the  stress 
distributions  in  the  four-ply  laminate  specimen  between 
results  from  a  uniform  160-element  model  shown  in  Figure 
(45)  and  those  from  a  locally  refined  288-element  model.  The 
288-element,  Figure  (4),  model  was  obtained  by  dividing  the 
two  elements  closest  to  the  free  edge  in  the  160-element 
model  into  10  elements  along  the  Z-direction.  Refinement  for 
both  the  cross-ply  and  the  angle-ply  specimens  were  studied. 

Figure  (46)  shows  the  distribution  of  R-stress  along  R5 
near  the  interface  of  0/90  layers.  It  is  clear  that  the 
singular  stress  behavior  near  the  free-edge  was  reproduced 
more  closely  by  the  288-element  model  and  did  not  show  well 
in  the  results  based  on  the  160-element  mesh.  Therefore,  a 
finite  element  mesh  refined  in  regions  of  steeply  varying 
stresses  is  necessary. 

Further  refinement  in  the  Z-direction  was  considered  in 
order  to  study  the  effect  of  mesh  refinement  on  some  stress 
distributions  where  poor  results  were  observed.  A  more 
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Figure  45:  160  Element  Model  (16x10) 
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refined  mesh,  416  elements  (16x26)  for  half  section 
analysis,  obtained  by  further  dividing  the  two  elements 
closest  to  the  free  edge  in  the  288-element  model  (16x18) 
into  10  elements,  was  used.  Results  from  the  axisymmetric 
analysis  were  compared  with  the  Q23  element  solutions 
obtained  using  the  original  mesh  size  which  corresponds  to 
the  (16x18)  elements  model.  Figure  (47)  indicates  the 
improvement  in  the  distribution  of  Rz-stress  along  R5 
specially  near  the  free  edge  where  the  traction-free 
boundary  condition  was  better  reproduced  by  the  finer  mesh. 
The  distribution  of  Z-stress  along  R4  is  shown  in  Figure 

(48)  .  Refinement  did  not  make  any  noticeable  difference  to 
the  distribution.  The  only  improvement  was  that  values  of 
stresses  closer  to  the  free  edge  were  available.  Figures 

(49)  ,  (50)  show  the  distribution  of  R-stress  at  R7  and  R8. 
Refinement  of  the  mesh  near  the  free  edge  did  not  have  much 
influence  on  the  calculated  values  of  the  stress  components. 
The  oscillatory  variation  in  the  stress  component  near  the 
free  edge  is  present  in  all  the  analyses. 

The  effect  of  mesh  refinement  on  the  distribution  of 
Tz,  Rz  and  Z-stresses  along  R4  of  [45/-45]s  laminates  is 
shown  in  Figures  (51)  ,  (52)  and  (53)  respectively.  The 
refined  mesh  near  the  free  edge  was  necessary  to  represent 
the  steeply  varying  stresses  in  this  vicinity  (i.e.,  near 
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Figure  48:  Effect  of  Mesh  Refinement  on  Distribution 
of  Z-Stress  Along  R4  of  [0/90]s  Laminate 
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Figure  50:  Effect  of  Mesh  Refinement  on  Distribution 
of  R-Stress  Along  R8  of  [0/90]s  Laminate 
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the  interface  of  45/-45  layers  and  near  free  edge)  .  The 
figures  indicate  that  the  stress  values  obtained  using  the 
refined  mesh  approached  the  traction-free  boundary 
condition. 

Figure  (54)  shows  the  effect  of  mesh  refinement  near 
the  free  edge  of  TR-stress  along  R8 .  The  stresses  varied 
steeply  near  the  corner  of  the  laminate  and  the  singular 
stress  behavior  was  reproduced  somewhat  better  by  the 
refined  mesh. 

The  effect  of  refinement  in  the  r-direction  was  also 
studied.  Figures  (55) ,  (56)  show  the  distribution  of  R- 
stress  along  R7  and  R8  of  the  cross-ply  laminate  for  N=48 
(each  physical  layer  was  subdivided  into  twelve  uniform 
sublayers) ,  and  they  indicate  that  the  result  of  the  refined 
mesh,  N=48  gave  better  results  then  those  obtained  for  N=16. 
The  values  at  R7  were  nearly  coincident  with  the  results  of 
the  Q23  analysis.  At  R8  there  was  considerable  improvement 
but  an  oscillation  of  stresses  existed  near  Z/B  =  0.6  to 
0.8. 


The  effect  of  mesh  refinement  in  the  r-direction  on 
through-thickness  distribution  of  stresses  near  the  free 
edge  for  N=32,  (each  physical  layer  was  subdivided  into  8 
sublayers)  and  N=16  is  described  in  Figures  (57)  through 
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Figure  55:  Effect  of  Mesh  Refinement  in  F-Direction 
on  Distribution  of  R-Stress  Along  R7 
of  [ 0/90 ] s  Laminate. 
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(61) .  The  continuous-traction  analysis  gave  values  of  the 
stress  field  at  the  interface  of  the  layers.  For  the 
axisymmetric  analysis,  an  estimate  of  the  stress  field  at 
the  interface  of  lamina  was  made  by  interpolating  from  the 
stress  field  computed  at  the  centers  of  the  sublayers.  Cubic 
spline  functions  were  used. 

Figures  (57),  (58)  show  through-thickness  distribution 
of  R  and  Rz-stresses  at  Z/B=0.99  of  cross-ply  laminate.  The 
steeply  varying  stresses  near  the  interface  of  the  0/90 
laminate  and  the  stress  free  condition  at  the  free  surface 
were  more  closely  reproduced  by  the  refined  mesh.  The 
percentage  difference  in  the  values  of  the  stresses,  at  the 
interface  of  0/90  laminate,  using  the  present  analysis  and 
using  the  Q23  element  were  10%  and  46%  for  R-stress  and  Rz- 
stress  respectively. 

Figures  (59) ,  (60)  and  (61)  show  through-thickness 
distribution  of  TR,  R  and  Rz-stresses  near  the  free  edge  of 
[45/-45]s  laminate.  The  figures  indicate  a  difference  in 
results  between  the  axisymmetric  analysis  and  the  Q23 
element  solution  near  interface  of  45/-45  layers.  However, 
the  refined  mesh  gave  satisfactory  results  near  the  free 
surface  where  a  traction-free  condition  exists.  At  the 
interface  of  the  45/-45  layers,  the  interpolated  values 
obtained  using  spline  fitting  of  the  axisymmetric  analysis 
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Figure  57:  Effect  of  Mesh  Refinement  in  R-Direction 
on  Through-Thickness  Distribution  of 
R-Stress  at  Z/B  =  0.99  of  [0/90]s  Laminate 
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Figure  58:  Effect  of  Mesh  Refinement  in  R-Direction 
on  Through-Thickness  Distribution  of 
Rz-Stress  at  Z/B  =  0.99  of  [0/90] 
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Figure  59:  Effect  of  Mesh  Refinement  in  R-Direction 
on  Through-Thickness  Distribution  of 
TR-Stress  at  Z/B  =  0.99  of  [45/-45] 
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Figure  60:  Effect  of  Mesh  Refinement  in  R-Direction 
on  Through-Thickness  Distribution  of 
R-Stress  at  Z/B  =  0.99  of  [45/-45] 
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Figure  61:  Effect  of  Mesh  Refinement  in  R-Direction 
on  Through-Thickness  Distribution  of 
Rz-Stress  at  Z/B  =  0.99  of  [45/-45Js 
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results  differed  by  16%,  26%,  10%,  as  compared  to 
approximate  stress  value  obtained  using  the  Q23  analysis, 
for  TR,  R,  and  Rz  stress  components  respectively. 

Another  interpolating  scheme,  used  in  order  to  estimate 
the  stress  field  at  the  interface  in  Figures  (57)  through 
(61)  was  studied.  Cubic  spline  functions  were  used  except 
for  the  two  stress  values  near  the  interface  where  quadratic 
functions  were  used.  The  two  quadratic  curves  were  assumed 
to  have  equal  and  opposite  curvature.  In  this  case,  the 
interpolating  polynomials  had  continuous  first  and  second 
derivatives  except  at  the  interface  where  there  is  a 
discontinuity  in  stress  gradients.  Figures  (62) ,  (63)  show 
through-thickness  distribution  of  R  and  Rz-stresses  at 
Z/B=0 . 99  for  cross-ply  laminate.  Using  the  interpolating 
scheme  described  above,  the  percentage  difference,  at  the 
interface  of  0/90  layers,  of  R-stress  values  for  the 
axisymmetric  analysis  compared  with  that  from  the  Q23 
element  analysis  increased  from  10%  to  16%  while  that  of  Rz- 
stress  decreased  from  46%  to  28%,  as  compared  to  the  spline 
fit  interpolation.  Figures  (64),  (65),  (66),  show, 
respectively,  through-thickness  distribution  of  Tr,  R,  Rz- 
stresses  near  the  free  edge.  The  percent  error,  at  the 
interface  of  45/-45  layers  decreased  from  16%  to  2%  for  Tr- 
stress  and  from  10%  to  9%  for  Rz-stress  but  increased  from 
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Figure  62:  Interpolation  of  Stresses  at  Interface:  R- 
Stress  Through-Thickness  distribution  at 
Z/B  =  0.99  of  [ 0/90 ] s  (Combination  of 
Cubic  and  Quadratic  Interpolating  Function) 
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Figure  63:  Interpolation  of  Stresses  at  Interface: 

RZ-Stress  Through-Thickness  distribution 
at  Z/B  =  0.99  of  [0/90]s  (Combination  of 
Cubic  and  Quadratic  Interpolating  Function) 


SH 


FREE  5URFPCE 


TR  STRESS  x  10E-06  (PSI) 

.00  0,140  0,80  1,20  1,60  2.00 


♦  REF  (11) 

*  32*18 

-  SPLINE  FIT  OF  DATA  INDICRTED  BY  it 
WITH  QURDRRT 1 C  FUNCTION  NEAR  INTERFACE 


Figure  64:  Interpolation  of  Stresses  at  Interface: 

TR-Stress  Through-Thickness  distribution 
at  Z/B  =  0.99  of  [45/-45]g  (Combination  of 
Cubic  and  Quadratic  Interpolating  Function) 
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Figure  65:  Interpolation  of  Stresses  at  Interface: 

R-Stress  Through-Thickness  distribution 
at  Z/B  =  0.99  of  [45/-45]s  (Combination  of 
Cubic  and  Quadratic  Interpolating  Function) 
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Figure  66:  Interpolation  of  Stresses  at  Interface: 

RZ-Stress  Through -Thickness  distribution 
at  Z/B  =  0.99  of  [ 4 5/ -4 6 ] 8  (Combination  of 
Cubic  and  Quadratic  Interpolating  Function) 
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26%  to  35%  for  R-stress,  as  compared  with  the  spline  fit 
interpolation . 

A  third  scheme  consisted  of  a  two-step  curve  fitting 
procedure.  The  stress  values  on  each  side  of  the  interface 
were  obtained  by  extrapolating  from  a  spline  fit  to  all  the 
points  other  than  the  interface.  The  two  curves  had  zero 
curvature  at  the  interface  and  gave  different  values  of  the 
stress  at  the  interface.  The  actual  stress  at  the  interface 
was  then  assumed  to  be  the  average  of  the  two  interpolated 
values.  Having  found  the  stress  at  the  interface,  a  spline 
fit  of  data  was  carried  out.  The  results  showed  that  there 
was  not  much  difference  in  the  interpolated  stress  value  as 
compared  to  the  previous  interpolation  scheme.  Figures  (67)  , 
(68)  show  the  distribution  of  Rz-stress  for  the  cross-ply 
laminate  and  TR-stress  for  the  angle  ply  laminate, 
respectively,  for  this  scheme. 
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Figure  67:  Interpolation  of  Stresses  at  Interface: 

RZ-Stress  Through-Thickness  distribution 
at  Z/B  =  0.99  of  [ 0/90] s  (Two-Step  Curve 
Fitting) 
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Figure  68:  Interpolation  of  Stresses  at  Interface: 

TR-Stress  Through-Thickness  distribution 
at  Z/B  =  0.99  of  [45/-45]s  (Two-Step  Curve 
Fitting) 
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4.4  TWENTY-TWO  ANGLE  PLY  LAMINATE 

4.4.1  Geometry.  Material  Properties  and  Finite  Element 
Mesh 

A  ring  specimen  with  fiber  orientation  [(25.5/- 
25.5) g/90] s  was  used  in  the  present  investigation.  The 
laminate  width  and  ply  thickness  were  1.0  in  and  0.00505  in 
respectively.  The  material  properties  were  [11]: 

E11  =  19-26  *  106  Psi 

E22  =  E3 3  =  1-32  *  106  psi 

G12  =  G13  =  G23  =  0.83  *  106  psi 

^12  =  ^13  =  ^23  =  0 ‘ 35 

Each  ply  was  modeled  by  a  single  element  through  its 
thickness  (N  =  22)  .  The  mesh  used  was  more  refined  near  the 
free-edge.  Full  section  was  used  for  the  axisymmetric 
analysis,  corresponding  to  308  (22x14)  mesh  as  shown  in 

Figure  (69) ,  for  half-section  discretization  of  the  ring. 
This  corresponds  to  154-element  model  for  a  quadrant  as  used 
with  the  continuous  traction  Q23  element  in  [11].  The  value 
R^,  shown  in  Figure  (69),  denotes  the  distance  from  the 
central  plane  of  section  to  the  center  of  the  ifcil  layer.  R3 
is  the  closest  to  the  midplane  of  the  specimen  and 
corresponds  to  the  90-degree  ply,  R-^  is  closest  to  the  top 
surface  and  corresponds  to  the  25.5-degree  ply.  The  internal 
pressure  applied  was  0.00987687  (psi)  which  yields  an 
average  axial  strain  equal  to  that  applied  in  Q23  element 
analysis  (0.95414*10“5) . 


4.4.2  Numerical  Results 


Numerical  results  based  on  the  pseudo-three  dimensional 
axisymmetric  analysis  using  QM4  element  were  compared  with 
results  based  on  the  continuous  traction  Q23  element 
described  in  [11].  The  same  mesh  size  was  used  in  both  the 
analysis.  Figures  (70)  through  (85)  show  comparison  of 
stress  fields  at  specific  locations  for  the  multi-ply 
laminate  using  the  Q23  element  and  the  axisymmetric 
analysis. 

Figures  (70),  (71),  (72)  and  (73)  show  the  distribution 
of  z-stress  along  Rl,  R5,  R6  and  Rll  respectively. 
Satisfactory  agreement  was  observed  between  the  axisymmetric 
analysis  and  the  Q23  element  solution  except  at  Z/B  =  0.56 
to  0.89  in  Figures  (70)  and  (73).  This  could  be  due  to  the 
finite  element  mesh  not  being  fine  enough  in  this  region. 
The  effect  of  mesh  refinement  is  described  later  in  this 
section.  The  figures  also  indicate  that  the  stress  values 
obtained  using  the  axisymmetric  analysis  approach  the 
traction-free  boundary  condition  at  the  free  edge. 

A  comparison  of  Rz-stress  distribution  along  Rl,  R5  and 
R6  (Figures  (74)  ,  (75)  and  (76)  respectively)  indicated  that 
the  results  obtained  using  axisymmetric  analysis  agreed  well 
with  Q23  element  solution  across  the  entire  width  of  the 
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Figure  70:  Distribution  of  Z-Stress  Along  R1  of 
[ (25.5/-25.5)s/90]s  Laminate 
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Figure  71:  Distribution  of  Z-Stress  Along  R5  of 
[  (25.5/-25.5)s/90]s  Laminate 
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Figure  72:  Distribution  of  Z-Stress  Along  R6  of 
[  (25.5/-25.5)s/90]s  Laminate 
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Figure  75:  Distribution  of  RZ-Stress  Along  R5  of 
[  (25.5/-25.5)s/90]s  Laminate 
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Figure  76:  Distribution  of  RZ-Stress  Along  R6  of 
[ (25.5/-25.5)s/90]8  Laminate 
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laminate  except  in  the  region  of  Z/B  =  0.56  to  0.79  where 
apparently  the  mesh  was  not  fine  enough  to  approximate  the 
steeply  varying  stress.  Figure  (77)  shows  the  distribution 
of  Rz-stress  at  Rll,  center  of  the  layer  nearest  to  the 
surface.  A  discrepancy  in  results  was  observed. 

Figures  (78),  (79),  (80)  and  (81)  show  the  distribution 
of  Tz-stress  at  Rl,  R5,  R6  and  Rll  respectively.  A  close 
agreement  was  observed  between  the  Q23  element  and  the 
axisymmetric  solution  except  near  the  free  edge  where  a 
steeply  varying  stress  exists.  Distribution  of  T-stress 
along  Rl  is  shown  in  Figure  (82) .  The  figure  indicates 
satisfactory  agreement  between  the  axisymmetric  analysis  and 
the  Q23  element  solution. 

Figures  (83) , (84) ,  (85)  show,  respectively,  through¬ 
thickness  distribution  of  R,  Rz  and  TR-stresses  at 
Z/B=0 . 995 .  A  difference  in  results  at  some  locations  was 
observed  between  the  two  analyses  for  the  Rz  and  TR- 
stresses.  The  R-stress  values  at  centers  of  elements  were 
underestimated  throughout.  Considering  that  the  coarseness 
of  mesh  could  be  contributing  to  the  discrepancy  and  to 
study  convergence  of  results  effect  of  mesh  refinement  in 
both  the  R  and  the  Z-direction  was  considered. 


Figure  77:  Distribution  of  RZ-Stress  Along  Rll  of 
[  (25. 5/-25. 5)S/90]S  Laminate 

118 


TZ  STRESS  x  1 OE-02  (PSD 

.05  0.13  0.01  0.29  0.37  0. 45 


REF  (11) 


o  PRESENT  ANALYSIS 


DISTANCE  FROM  CENTER  LINE  Z/B 


Figure  81:  Distribution  of  TZ-Stress  Along  Rll  of 
[(25.5/-25.5)5/90]#  Laminate 
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Figure  82:  Distribution  of  T-Stress  Along  R1  of 
[  (25.5/-25.5)s/90]s  Laminate 
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Figure  83:  Through-Thickness  Distribution  of  R 
Stress  Near  Free  edge  (Z/B=0.995) 
of  [  (25.5/-25.5)5/90]g  Laminate 
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Figure  84:  Through-Thickness  Distribution  of  Rz- 
Stress  Near  Free  edge  (Z/B=0.995) 
of  [ (25.5/-25.5)s/90]s  Laminate 
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Figure  85:  Through-Thickness  Distribution  of  TR 
Stress  Near  Free  edge  (Z/B=0.995) 
of  [(25.57-25.5)5/90]  Laminate 


4.4.3  Effect  of  Mesh  Refinement 


Mesh  refinement  in  the  R-direction  as  well  as  in  the  Z- 
direction  near  the  free  edge  was  considered.  Full  section 
was  used  for  the  axisymmetric  analysis.  The  mesh, 
corresponding  to  726  (22x33)  as  shown  in  Figure  (86)  for 
half  section  discretization,  was  used  to  study  refinement 
effect  near  the  free  edge.  The  mesh  was  obtained  by  dividing 
the  two  elements  closest  to  the  free  edge  in  the  308-element 
model  into  10  elements  and  also  by  refining  in  the  region  of 
Z/B  =  0.4  to  0.8.  results  from  the  axisymmetric  analysis 
were  compared  with  the  Q23  solutions  obtained  using  the 
original  mesh  size  which  corresponds  to  the  308-element 
model . 


Figures 

(87),  (88), 

(89)  ,  and  (90)  show 

the 

distribution 

of  Z-stress 

along  Rl, 

R5 ,  R6 ,  and 

Rll 

respectively . 

The  figures 

indicate 

improvement  in 

the 

distribution  of  Z-stress  specially  in  the  region  of  Z/B  = 
0.56  to  0.8  and  near  the  free  edge  where  the  traction-free 
boundary  condition  was  better  reproduced  by  the  finer  mesh. 

The  effect  of  mesh  refinement  on  Rz-stress  distribution 
along  Rl,  R5,  R6  and  Rll  is  shown  in  Figures  (91)  ,  (92)  , 
(93)  and  (94)  respectively.  The  refined  mesh  in  the  region 
of  Z/B  =  0.6  to  0.8  was  necessary  to  represent  the  stress 
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Figure  87:  Effect  of  Mesh  Refinement  on  Distribution 
of  Z-Stress  Along  Ri  of  [  (25/-25) ,/90] 
Laminate  5  s 
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Figure  88:  Effect  of  Mesh  Refinement  on  Distribution 
of  Z-Stress  Along  R5  of  [ (25/-25)  =/90] 
Laminate  8 
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Figure  89:  Effect  of  Mesh  Refinement  on  Distribution 
of  Z-Stress  Along  R6  of  [ (25/-25) /90] 
Laminate  8 
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Figure  90:  Effect  of  Mesh  Refinement  on  Distribution 
of  Z-Stress  Along  Rll  of  [ (25/-25) j/90],, 
Laminate 


132 


FREE  EDGE 


-  REF  (11) 
O  22x14 
A  22x33 


Figure  91:  Effect  of  Mesh  Refinement  on  Distribution 
of  RZ-Stress  Along  R1  of  [ (25/-25) c/90] 
Laminate  '  * 
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variation  in  this  region  even  though  a  difference  in  the 
axisymmetric  analysis  and  the  Q23  element  solution  was  still 
observed  in  Figures  (92-94) .  The  distribution  of  Rz-stress 
along  Rll  (Figure  (94))  shows  that  the  stresses  varies 
steeply  near  the  corner  of  the  laminate  and  the  singular 
stress  behavior  was  reproduced  somewhat  better  by  the 
refined  mesh. 

Figures  (95)  ,  (96)  ,  (97)  show  the  effect  of  mesh 
refinement  on  Tz-stress  distribution  along  Rl,  R5,  and  R6. 
A  difference  in  results  between  the  Q23  element  solution  and 
the  axisymmetric  analysis  near  the  free  edge  was  still 
observed.  Figures  (98)  shows  the  distribution  of  Tz-stress 
along  Rll  and  indicates  that  the  results  obtained  using  the 
refined  mesh  agreed  well  with  the  Q23  element  solution. 

The  effect  of  mesh  refinement  in  the  R-direction  on 
through-thickness  distribution  of  stresses  near  the  free 
edge  was  also  studied.  Each  physical  layer  was  subdivided 
into  two  sublayers,  which  correspond  to  (44x14)  element 
model  for  half  section  discretization.  For  the  axisymmetric 
analysis,  an  estimate  of  the  stress  field  at  the  interfaces 
of  the  lamina  was  made  by  interpolating  from  the  stress 
field  computed  at  the  adjacent  centers  of  the  sublayers. 
Cubic  spline  functions  were  used.  Results  were  compared  with 
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Figure  96:  Effect  of  Mesh  Refinement  on  Distribution 
of  TZ-Stress  Along  R5  of  [ (25/-25) 5/90]s 
Laminate 
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Figure  97:  Effect  of  Mesh  Refinement  on  Distribution 
of  TZ-Stress  Along  R6  of  [ (25/“25) 5/90] s 
Laminate 
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the  Q23  element  solution  which  gave  values  of  the  stress 
fields  at  the  interfaces.  Figure  (99)  shows  through¬ 
thickness  distribution  of  R-stress  at  Z/B=  0.995.  The  Q23 
element  solution  showed  an  oscillatory  pattern  of  stress 
which  was  not  seen  in  the  spline  fit  of  the  axisymmetric 
analysis.  Figure  (100)  shows  through- thickness  distribution 
of  Rz-stress  at  Z/B  =  0.995.  The  figure  indicates 
satisfactory  agreement  between  the  axisymmetric  analysis  and 
the  Q23  element  solution  except  at  the  interface  of  the  90- 
degree  layer  and  the  angle-ply  layer  where  a  singular  stress 
behavior  is  observed.  Figure  (101)  shows  through-thickness 
distribution  of  Tr-stress  near  the  free  edge.  The  refined 
mesh  showed  an  oscillatory  pattern  of  stress  which  agreed 
with  the  Q23  solution.  At  the  interfaces  of  the  layers,  the 
interpolated  value,  obtained  using  spline  fitting  of  the 
axisymmetric  analysis  results  differed  by  up  to  16%  as 
compared  to  the  approximate  stress  values  obtained  using  the 
Q23  analysis.  However,  this  discrepancy  could  be  largely  due 
to  the  spline  fit  which  insures  continuity  of  slopes  and 
curvatures.  For  the  laminate  considered,  the  TR-stress 
component  may  (and  Q23  analysis  confirms  this)  have 
discontinuity  in  gradients  at  the  interfaces. 
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Figure  99:  Effect  of  Mesh  Refinement  in  R-Direction 
on  Through-Thickness  Distribution  of 
R-Stress  at  Z/B  =  0.995  Of  [ (25/-25) 5/90] 
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Figure  101:  Effect  of  Mesh  Refinement  in  R-Direction 
on  Through-Thickness  Distribution  of  TR- 
Stress  at  Z/B  =  0.995  of  [ (25/-25) 5/90]_ 


SECTION  V 


DISCUSSION 


The  problem  of  linear  elastic  stress  analysis  of  a 
composite  structure  of  ax i symmetric  geometry  and  loading  has 
been  investigated  and  computer  programs  developed  based  on 
finite  element  discretization  using  the  modified  bilinear 
quadrilateral  elements  with  selectively  reduced  integration. 
Allowing  for  the  tangential  degree  of  freedom,  the  analysis 
is  applicable  to  angle  ply  laminates.  The  axisymmetric 
analysis  procedure  can  be  used  to  model  behavior  of  cross- 
ply  as  well  as  angle-ply  coupons  subjected  to  uniform  axial 
strain.  For  a  sufficiently  large  diameter,  a  ring  specimen 
subjected  to  internal  pressure  is  equivalent  to  an  endless 
straight  free  edge  coupon  subjected  to  an  essentially 
uniform  inplane  extension.  Therefore,  the  main  advantage  of 
the  ring  specimen,  compared  with  other  finite  element 
models,  is  that  it  is  completely  free  of  end  constraints, 
thereby  modeling  the  free  edge  delamination  coupons  somewhat 
better. 
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In  the  cross-ply  laminate  analysis,  numerical  results 
calculated  from  the  axisymmetric  stress  analysis  generally 
agreed  well  with  the  continuous  traction  Q23  element  except 
near  the  free  surface.  This  is  because  the  continuous 
traction  Q23  element  allows  traction-free  boundary 
conditions  to  be  specified  in  a  pointwise  sense  while  the 
present  analysis  does  not.  The  study  also  revealed  that  the 
pattern  of  mesh  refinement  had  significant  effect  on  the 
estimate  of  the  stress  distribution  in  the  vicinity  of  the 
free  surface  and  the  traction-free  edge. 

For  the  angle-ply  laminate  problem,  a  pseudo-three 
dimensional  analysis  specifying  geometry  and  displacements 
as  axisymmetric  but  including  displacements  in  all  three 
coordinates  directions  is  necessary.  Further,  discretization 
of  the  full  section  is  necessary  for  the  analysis  of  angle- 
ply  laminate  since  symmetry  about  the  Z  =  0  plane  does  not 
exist.  Numerical  results  from  the  axisymmetric  analysis 
agreed  well  with  continuous  traction  Q23  element  except  near 
the  free  edge.  Mesh  refinement  in  the  vicinity  of  the  free 
edge  indicated  that  the  traction-free  boundary  condition  was 
more  closely  approximated  by  the  refined  mesh.  It  was  found 
that  for  a  good  estimate  of  the  stress  field,  a  mesh  refined 
in  the  vicinity  of  the  free  edge  as  well  as  near  the  free 
surface  should  be  used. 


In  the  multi-ply  laminate  analysis,  numerical  results 
obtained  using  the  pseudo-three  dimensional  axisymmetric 
analysis  agreed  well  with  the  Q23  element  solution  except  at 
the  free  edge  and  at  the  interfaces  of  the  layers.  A  mesh 
refined  in  the  vicinity  of  the  free-edge  and  in  the 
thickness  direction  had  significant  effect  on  the  estimate 
of  the  stress  field. 

The  computational  procedure  developed  is  efficient  in 
terms  of  computation  time  for  the  analysis  of  axisymmetric 
composite  structures  with  axisymmetric  loading.  In 
particular,  for  the  free-edge  delamination  problem,  the 
procedure  is  much  less  expensive  than  the  continuous 
traction  Q23  element  [11].  For  the  same  mesh  of  the  four-ply 
laminate  problem,  the  CPU  time,  on  an  IBM  3081  mainframe 
computer,  using  the  present  analysis  was  25  seconds  for  the 
two-dimensional  analysis  (323  nodal  points),  and  4  minutes 
for  the  pseudo-three-dimensional  analysis  (629  nodal  points) 
and  full  section  discretization,  as  compared  to  40  minutes 
using  the  continuous  traction  Q23  element  (485  nodal 
points)  .  For  the  22-ply  laminate  problem,  the  CPU  time  using 
the  present  analysis  was,  respectively,  4  minutes  40 
seconds,  9  minutes  41  seconds,  12  minutes  and  31  minutes  45 
seconds  for  a  mesh  corresponding  to  (22x14)  (667  nodal 
points),  (44x14)  (1305  nodal  points),  (22x33)  (1541  nodal 
points)  and  (44x33)  (3015  nodal  points)  for  half  section 
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discretization.  The  CPU  time  using  the  continuous  traction 
element  was  45  minutes  38  seconds  for  a  mesh  of  (14x22)  (513 
nodal  points)  .  The  code  is  applicable  to  general 
axisymmetric  laminated  systems. 
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